Traveling Waves, Front Selection, and Exact Nontrivial Exponents in a Random 

Fragmentation Problem 
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We study a random bisection problem where an initial interval of length x is cut into two random 
fragments at the first stage, then each of these two fragments is cut further, etc. We compute the 
probability P n (x) that at the n th stage, each of 2 n fragments is shorter than 1. We show that P n (x) 
approaches a traveling wave form, and the front position x n increases as x n ~ n 13 p" for large n. 
We compute exactly the exponents p = 1.261076 . . . and /3 = 0.453025 ... as roots of transcendental 
equations. We also solve the m— section problem where each interval is broken into m fragments. 
In particular, the generalized exponents grow as p m ~ m/(lnm) and /3 m ~ 3/(2 In m) in the large 
m limit. Our approach establishes an intriguing connection between extreme value statistics and 
traveling wave propagation in the context of the fragmentation problem. 



PACS numbers: 02.50.-r, 05.40.-a, 64.60.-i 

The statistics of extremes plays an important role in 
various branches of physics, statistics, and mathemat- 
ics For example, in physics of disordered sys- 
tems, the statistics of extremely low energy states gov- 
erns the thermodynamic behavior at low temperatures 
H . The extreme- value statistics is well understood when 
random variables are independent and identically dis- 
tributed. Much less is known when the random variables 
are correlated. In the replica language, this class of prob- 
lems corresponds to full replica symmetry breaking 
Thus it would be important to derive exact results for 
extreme value statistics for correlated random variables. 

On the other hand, there is a wide variety of problems 
in physics, chemistry and biology, that allows for solu- 
tions with propagating traveling waves Q| . The problem 
of front propagation into an unstable state dates back 
to the pioneering studies motivated by gene spread- 
ing in a population. Numerous traveling wave solutions 
were found in combustion |7j and reaction diffusion sys- 
tems ||; they also appear in the mean-field theory of 
directed polymers in random medium |^] , calculations of 
Lyapunov exponents of random matrices Jlfj|,pd| , pattern 
formation and many other problems ]13[|. Usually 
in all these problems, the traveling wave has an expo- 
nentially decaying front that advances with a uniform 
velocity from a stable state to an unstable state as time 
increases. Out of a continuum of possible allowed veloci- 
ties of the front, a unique value is usually selected. This 
velocity selection mechanism has been investigated by a 
variety of methods p.pp^ |i"6| ] , and it was found that for 
"sufficiently steep" initial conditions, the minimal veloc- 
ity is usually selected. 

Is there a connection between these two sets of prob- 
lems? In this paper, we give a positive answer to this 
question in the context of a random fragmentation prob- 
lem jl7| . Our study indicates that there is perhaps a very 



general and deep connection between the two seemingly 
unrelated topics of extreme value statistics and propagat- 
ing wave solutions. Besides, we give an exactly solvable 
example of extreme value statistics for correlated random 
variables. 

The random bisection problem can be formulated as 
follows. An interval of length x is cut into two halves of 
lengths x\ = rx and X2 = (1— t)x, respectively, where r is 
the random number chosen from the uniform distribution 
over [0,1]. Next, each of these two fragments is cut again 
cut into two parts (with no correlations between the cut- 
tings). After the n th step, there are 2™ fragments whose 
lengths are correlated random variables (correlations are 
dynamically generated). Given the initial length x, what 
is the probability P n (x,l) that each of the 2™ fragments 
is shorter than 11 In other words, what is the probability 
that the longest of these 2™ fragments would be shorter 
than I. The I dependence is simple, P n (x,l) = P n (x/l), 
so in the following we set I = 1 without loss of generality. 

Clearly, P n (x) = for x > 2™. Moreover, there ex- 
ists a threshold x n such that for large n, P n (x) — > 1 for 
x <C x n and P n (x) — > for x 3> x n . The structure is 
that of an advancing front and the front position x n in- 
creases with n as p n for large n. This intriguing "phase 
transition" was noticed by several workers |l8| and 
finally rigorously proved by Devroye pH who also com- 
puted exactly the value of p — 1.261076.... Recently, 
Hattori and Ochiai p2| conjectured, based on numerical 
simulation, that the more precise asymptotic behavior of 
x n reads 



(1) 

and found numerically that (3 « 0.407. 

In this paper, we derive this asymptotic result by using 
methods of statistical physics. In particular, we re-derive 
the value of p in a physically transparent way and com- 
pute the exponent [3 exactly, which is a completely new 
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result. We show that (3 = 0.453025 ... is nontrivial and 
is given by the root of a transcendental equation. Our 
technique consists of employing a scaling ansatz which 
reduces the problem to solving an equation that admits 
traveling wave solutions. We then select an appropriate 
solution by employing a well known front selection princi- 
ple. Our technique also allows us to obtain exact results 
for a multisection problem where each interval is cut into 
to random pieces at every stage. Once again, the distri- 
bution P n (x) follows the zero-one law, with the threshold 
value x n (m) ~ n^ m [p m ] n . As in the to — 2 case, we com- 
pute the exponents p m and (3 m exactly for arbitrary to. 
For example, p 3 = 1.499118... and = 0.429815..., 
and p m — ► m/lnm and f3 m — > 3/(2 In to) when to — > oo. 
We also show that the variance of the front position 
x n (m) is always finite. 

We first consider the m = 2 case, i.e., the bisection 
problem. It is easy to see that P n (x) satisfies the exact 
recursion relation, 



P 



n+l{x) = ~ 
X 



dyP n (y)P n {x - y), 



(2) 



with the initial condition, Pq(x) = 9(1 — x), where 9(x) 
is the Heaviside step function. The prefactor 1/x on the 
right hand side of Eq. (||) is the probability that the point 
where the interval is cut into two pieces is chosen ran- 
domly from the interval [0, a;]. Taking the Laplace trans- 
form Q n (s) = J °° dxe~ sx P n (x) of Eq. (|), we get, 



dQn+l(s) 

ds 



-Ql(s). 



(3) 



Physically, one expects that as n grows, P n (x) will be 
nonzero for x < x n and then will rapidly decay to zero 
for x > x n where threshold. An appropriate def- 

inition of x„ would be, 



dx P n (x). 



(4) 



Thus, it is natural to rewrite P n (x) in the scaling form, 
Pn[x) = f n (x/x n ). Thence, the Laplace transform reads 
Q n (s) = x n F n (sx n ) with F n (z) = J^°jlye 
Substituting this form of Q n (s) into Eq. 



dF n+1 (z) 
dz 



X n +1 



2 r 



F„ 



X n +1 



zy fn(y)- 
we get 

2 

(5) 



Note that, by definition, Q n (0) = x„ implying F n (0) = 1 
for all n. Additionally, / n (0) = 1 leads to F n (0) — > z' 1 as 
z —> oo. For convenience, we make a further substitution, 
F n (z) = (1 — H n (z))/z, which recasts Eq. (0) into 



dH n 



+i 



dz 



X n +1 



2H„ 



X n +1 



H n (z), (6) 



with the boundary conditions H n (0) 
as z — > oo. 



1 and H n (z) -> 



We first focus on the computing of p in the asymptotic 
relation (|l]). To determine p, we seek a 'stationary', i.e., 
independent of n, solution of Eq. (^|). Using x n ~ n^p 11 
and taking the n — > oo limit, we find that the stationary 
solution satisfies 



dz 



- 2H 



H{z), 



(7) 



where H(0) = 1 and H(z) — > as z — > oo. While we 
could not solve this non-local and nonlinear differential 
equation, we can determine the "eigenvalue" p exactly 
through the asymptotic analysis. Indeed, in the large z 
limit, H(z) is small and thus one can neglect the nonlin- 
ear term in Eq. (Q). The resulting linear equation admits 
a power law solution, H(z) = az 



A with 



P = 



(8) 



Thus, a wide range of possible /i's is in principle allowed. 
However, usually a particular value is selected depending 
on the initial condition of the system. This is very sim- 
ilar to the problem of velocity selection in a large class 
of problems with wave propagation [p|JT^] and it is well 
known that for a wide class of initial conditions, the ex- 
tremum value is chosen. In the present case, the function 
on the right hand side of Eq. (^) has a unique maximum 
at fx = (J,*, which is a root of ln( + £ ) = p*/(l + /-**)• 
Though we have not proved explicitly that the extremum 
value is indeed chosen, one can infer this conclusion from 
the general principle of front selection. Note also that 
p = p(p*) = 1.261070 . . . can be written as p — e a where 
a is a solution of a = log(2ae), in agreement with the 
result of Ref. fH|. 

The exponent [3 characterizes the next to leading 
asymptotic behavior of the front. Therefore, to compute 
(3 we need to consider the full equation (|^) rather than its 
n — > oo limit. A sub- leading asymptotic behavior of trav- 
eling fronts was originally analyzed by Bramson ]l4| for a 
reaction-diffusion equation, and recently investigated in 
]l5| , p~6 13 1 . Here we employ an approach of Ref. [^5). For 
finite but large n, we make the following scaling ansatz 
for the function H n (z), 



H n {z) 



n c 



(9) 



where the exponent a and the scaling function G(y) are 
yet to be determined. The scaling function G(y) must 
vanish as y — * oo. Also, G(y) ~ y as y — > to en- 
sure that for large n, H n (z) ~ and is independent 
of n. We substitute this scaling ansatz in Eq. (||) and 
use x n ~ rfi p n where p is already known exactly from 
Eq. (||). Using the exact value of (i*, we find that dif- 
ferent leading order terms are comparable only with the 
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special choice a = 1/2. In that case, the scaling function 
G(y) satisfies an ordinary differential equation, 



d 2 G 



dG 



{\npY — + y— + (2/3// - l)G(y) = 0, 



dy 



dy 



(10) 



with the boundary conditions G(y) ~ y for y — > and 
G(y) — > as y — > oo. This therefore constitutes an eigen- 
value problem where (3 is the required eigenvalue. The 
exact solution of this differential equation that satisfies 
the boundary condition at y — > oo is given by [E3J , 



G(y) = Ae 



D 



y 

In p 



(11) 



where A is a constant and D p {x) is the parabolic cylinder 
function of order p and argument x. From the known 
properties of cylinder functions |p3| , we find that the 
boundary condition, G(y) ~ y as y — > 0, selects the eigen- 
value 2(/3/i* — 1) = 1. This is because only the parabolic 
cylinder function with index 1 vanishes linearly with x 
as x — * 0. Thus, the exponent (3 is exactly determined in 
terms of the known /i*, 



2p* 



(12) 



Using p* = 3.311070 . . ., we get (3 = 0.453025 . . .. Our 
exact result slightly differs from the numerical value 
j3 ks 0.407 obtained in Ref. An accurate numeri- 

cal determination of (3 is not simple as it is a sub-leading 
correction to the the leading asymptotic behavior. 

The bisection problem can be straightforwardly gen- 
eralized to the m-section problem where at each stage, 
every interval is cut into m random pieces [p4| . The prob- 
ability P n (x) that at the n th stage each of m™ fragments 
is shorter than 1 satisfies the exact recursion relation, 

= f .../ Y[d yj p n ( yj ) 

X Jo Jo J = 1 



Slj^yi-x), (13) 



with the initial condition, Pq(x) = 9(1 — x). This equa- 
tion can be easily derived as follows. At a given stage 
n, an interval of length x is cut into m fragments. Let 
z\, . . . , z m -i denote the location of the points at which 
the interval is cut. The allowed range of values of the co- 
ordinates [zi, . . . , 2 m _i] is the (m — l)-dimensional sim- 
plex: x > z m -i > . . . > Z\ > 0. The volume of this 
simplex is x m ~ Y /(m — 1)!, and this explains the pref- 
actor of Eq. (13). Finally, by changing co-ordinates to 
Vi = zi,y 2 = z 2 - Zi, . . . , y m = x - z m _ l7 one obtains 
Eq. (^). Note also that for m = 2, Eq. ( ]l3| ) reduces to 
Eq. (3). We now proceed exactly as in the m = 2 case. 
The Laplace transform, Q n (s), satisfies the differential 
equation, 



L Qn+l(.s) 



ds 171 ' 1 



= (-ir- 1 (m-l)\[Q n (s)f 



(14) 



We again expect the threshold to grow as n l3m (p m ) n . To 
determine p m , we assume that P n approaches the scaling 
form, P n (x) — f(x/x n (m)), in the large n limit. Thence, 
Q n (s) = x n (m)F(sx n (m)) where x n (m) = dx P n (x) 
and F(z) = J °° dz e~ zy f(y). By inserting the scaling 
form for Q n (s) and x n (m) ~ n l3m (p m ) n into Eq. ( p^| ) we 
finally arrive at the non-local differential equation, 



- x F(z) _ (-l) m - x (m - 1)! 



dz' 



(PmY 



Pm 



(15) 



Substituting F = (1 — H(z))/z and linearizing the re- 
sulting equation for H(z) in the large z limit, one finds a 
solution that behaves algebraically, H(z) ~ z^^ m , in the 
large z limit. A straightforward algebra then shows that 
p m depends on p m via 



r(Mr, 



r( Mm + i)r(m + 1) 



(16) 



where T(x) is the usual Gamma function. Equation ( |l6| ) 
is the generalization of Eq. (^) for the arbitrary m case. 
Once again, the function on the right hand side has a 
unique maximum at p m = [i* m and this maximum is ac- 
tually selected. In particular, p 3 — exp[ 
b is found from 6(26+3) = (6+l)(6+2) 



(2b+3) 



where 



L (6+1) (6+2) J 

lln[(6+l)(6+2)/6]. 
Solving this numerically, gives ps = 1.499118... and 
p* 3 = 3.489870 . . .. 

The exponent (3 m can be determined exactly for arbi- 
trary m. We do not repeat the calculation as it follows 
the same steps as for m = 2 and the final result is given 
by the same expression (|T2]), i.e., (3 m — i/(2p,* m ). For 
instance, (3% = 0.429815 . . ., and generally the exponent 
(3 m decreases with increasing m. 

One can easily derive the asymptotic behavior of p m 
and P m for large m. Taking logarithms on both sides 
of Eq. ( |T6| ) and differentiating with respect to p, m we 
find, after straightforward algebra, that the maximum 
occurs at p* m « lnm. In this calculation, we have 
used properties of Gamma function for large arguments: 
T(m + p m )/T(m + 1) ~ (m + IV*™- 1 for m > p, m . Sub- 
stituting p m = lnm into Eqs. ( |l6|) and (12), we find that 
to leading order for large m, 



lnm 



and f3 n 



2 In m 



(17) 



In this work, we have established a relationship be- 
tween two seemingly disparate subjects - the statistics 
of extremes and traveling wave propagation. Specifically, 
we have shown that the probability density of the maxi- 
mal fragment length, which is the derivative of the distri- 
bution P n (x), approaches the solitary traveling wave in 



3 



the large n limit. This traveling wave has a finite width 
which implies, in the context of the random binary search 
trees, that the variation of the height of a tree is finite. 
This result has long been anticipated on the basis of nu- 
merical experiments but remained intractable |2^,^5| . 

The present work admits several extensions. One could 
assume that an interval can be cut into fragments of rela- 
tive lengths r and 1 — r with the probability density 7r(r) 
which is arbitrary apart from normalization and symme- 
try requirements, 7r(r) — n(l — r) and ^drn{r) = 1. 
Besides the uniform probability density, 7r(r) = 1, one 
could solve the random bisection problem for a number 
of other densities, e.g., for 7r(r) = 6r(l — r). The basic re- 
sult, Eq. (Q), always holds but the exponents do depend 
on the probability density 7r(r). 

One could also modify the random bisection problem 
by deciding to cut only those intervals which are still 
longer than 1 . This problem can be solved by employing 
similar techniques as the original one. The most inter- 
esting question is how does the total number of intervals 
N(x), which are left after all intervals will be shorter 
than 1, depend on the initial length x. For the orig- 
inal random bisection problem with the uniform prob- 
ability density, our solution implies N(x) ~ x 1 (lnx)~ s 
with 7 = In 2/ ln/j = 2.9881 and 6 = 07 = 1.35368. For 
the modified version, we have found N(x) ~ 2x. Thus, 
the modified algorithm is much more effective if we want 
to minimize the total number of cuts. 

Apart from the straightforward applications of our re- 
sults to the random search tree problem in computer sci- 
ence, our results have possible implications in a number 
of topics of current interest in physics and chemistry. One 
obvious application is to the fracture of a rock, or to the 
beak-up of a polymer. Another possible application is to 
granular materials. Recent experiments have studied the 
propagation of stresses in a granular pile of glass beads 
subjected to a large vertical overload Our model is 
closer to the situation when the vertical overload is local- 
ized. This force gets transmitted from the top layer to 
the bottom layer. If this external force greatly exceeds 
the weight of individual grains, F 3> w, and if fractions 
of force transmitted from a grain to its neighbours in the 
lower layer are random, then the model studied above 
gives the mean-field description of the force transmission. 
Within this approximation, the force chains starting from 
the top never intersect with each other thus maintaining 
the tree structure. Our results then imply that if the 
force from a grain in a given layer is always transmitted 
to m grains in the next layer, than the granular material 
should contain at least ln(F/F*)/ In p m layers to guaran- 
tee that the normal force supported by any grain at the 
bottom layer never exceeds . 
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